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We numerically investigate the relaxation dynamics in an isolated quantum system of interacting 
bosons trapped in a double-well potential after an integrability breaking quench. Using the statistics 
of the spectrum, we identify the postquench Hamiltonian as nonchaotic and close to integrability over 
a wide range of interaction parameters. We demonstrate that the system exhibits thermalization in 
the context of the eigenstate thermalization hypothesis (ETH). We also explore the possibility of 
an initial state to delocalize with respect to the eigenstates of the postquench Hamiltonian even for 
energies away from the middle of the spectrum. We observe distinct regimes of equilibration process 
depending on the initial energy. For low energies, the system rapidly relaxes in a single step to a 
thermal state. As the energy increases towards the middle of the spectrum, the relaxation dynamics 
exhibits prethermalization and the lifetime of the metastable states grows. Time evolution of the 
occupation numbers and the von Neumann entropy in the mode-partitioned system underpins the 
analyses of the relaxation dynamics. 
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I. INTRODUCTION 

Motivated by experiments in ultracold atoms [1, 2], 
there has been a resurgence of interest in exploring the 
fundamental aspects of statistical mechanics and whether 
statistical properties can emerge from unitary time evo¬ 
lution. To this end, recent attempts in describing ther¬ 
malization in closed quantum systems focused on time- 
evolving observables that are local to a subsystem and 
how they approach thermal equilibrium described by 
Gibbs ensemble [3, 4]. One proposed mechanism to un¬ 
derstand thermalization is the eigenstate thermalization 
hypothesis (ETH) [4-6]. 

The ETH conjectures that under certain conditions, 
the expectation value of a physically relevant observable 
will have a long-time average close to an appropriate 
microcanonical ensemble prediction. When testing the 
ETH, the system is usually driven out of equilibrium by 
a quantum quench in the parameters of the Hamiltonian. 
The ETH was numerically tested against other hypothe¬ 
ses and it was shown that for initial states following the 
conditions of the ETH, thermalization is achieved within 
observables [4, 7]. To introduce basic ideas of the ETH, 
consider an initial state |(/)o), which is not an eigenstate 
of some postquench or final Hamiltonian H without any 
degeneracies. Then, the initial state can be expanded in 
terms of the eigenstates of H with eigenvalues as 

k 

After the quench, the initial state will undergo time evo¬ 
lution as 

\m = ( 2 ) 


and therefore, the expectation value of an observable is 
given by 

(0(0) = ^ (3) 

k,l 

with Oik = (^|0|fc). Due to dephasing, Eq. (3) will relax 
to the infinite-time average, 

f dT{d{T)) = V \ak\^Okk- (4) 
t Jq ^ 

This long-time average is commonly referred to as the 
diagonal ensemble prediction due to the diagonal ensem¬ 
ble, pd = (^1- According to the ETH, the 

diagonal ensemble average of an observable, Eq. (3), is 
close to an appropriate microcanonical average (0)me, 

(0)me=A/'“^ ^ Okk, (5) 

l-E'fc —-E'o|<A£J 

k 

(0)me = (^, 

where Nf is the number of eigenstates within a narrow 
energy window AE around the initial mean energy of 
the system Eq- The last equality in Eq. (5) is true if the 
eigenstate expectation value (EEV), Okk, is a smooth 
function of Ek, while the off-diagonal elements Oki are 
negligible [4, 6]. There is no rigorous analytical proof of 
the ETH but different numerical studies on spin systems 
and lattice models both for small and large system sizes 
support its validity [8-16]. 

Another important aspect of relaxation is the ques¬ 
tion of how local observables equilibrate. This is par¬ 
ticularly important in the dynamics of nearly integrable 
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systems, which may exhibit a two-stage equilibration pro¬ 
cess. Both theoretical and experimental studies have 
investigated the nature of prethermalization in various 
physical setups [17-27]. 

The emergence of statistical relaxation can also be 
linked to the onset of quantum chaos as highlighted by 
previous studies [28-30]. In this context, particular em¬ 
phasis is laid on the form of the energy distribution of 
an initial state and how it compares with predictions of 
quantum chaos. We utilize this approach to demonstrate 
chaotic initial states even for energies far from the mid¬ 
dle of the spectrum of eigenenergies, an extension of the 
findings in Ref. [31]. 

In this work, we study the nonequilibrium dynamics 
of interacting bosons in a one-dimensional double-well 
potential, which can be modelled using a two-site mul¬ 
tilevel Bose-Hubbard Hamiltonian [32, 33]. Contrary to 
previous studies dealing with few atoms in multiple lat¬ 
tice sites, we study the other limit of having finite large 
number of bosons confined in few sites. As a matter of 
fact, it was checked that this system will thermalize due 
to the ergodicity of the associated mean-field trajecto¬ 
ries in the semiclassical limit of large N [34]. Here, we 
only consider TV = 35 bosons and use exact diagonal- 
ization to access the eigenvalues and the eigenstates of 
the final Hamiltonian. First, we determine the chaoticity 
of the postquench Hamiltonian on the basis of spectral 
statistics. We then verify the validity of the ETH in this 
system. Our results based on the distribution of EEV 
illustrate that not all eigenstates are thermal. Using all 
possible initial product state configurations, we compare 
the expectation values of local operators obtained using 
the diagonal and the microcanonical ensemble. 

Most of the studies related to subsystem thermaliza- 
tion involve spatial partitioning of the system [12, 35-37]. 
Instead, we implement a different partitioning scheme 
where a subsystem will consist of only one mode, either 
the lower or the upper level in one of the wells. Since ther- 
malization is a dynamical process, we also compute the 
time evolution of local operators and the von Neumann 
entropy of a subsystem. We find that the relaxation dy¬ 
namics of the von Neumann entropy is instructive in re¬ 
vealing prethermalization in the system. Moreover, we 
observe that fast relaxation of a local operator correlates 
with a rapid single-step growth of the von Neumann en¬ 
tropy as it approaches the Gibbs entropy. On the other 
hand, we observe a two-stage relaxation process for initial 
states close to the middle of the spectrum. Finally, we 
discuss a possible mechanism for the prethermalization 
in the system and the associated time scales. 

This paper is organized as follows. In Sec. H, we first 
describe the model and the quench protocol used in this 
work. Section HI contains results that characterize the 
spectrum and the chaoticity of the postquench Hamilto¬ 
nian. We investigate the requirements of the ETH in Sec. 
IV. The key results on relaxation and thermalization in 
the system are discussed in Sec. V. In Sec. V A, we com¬ 
pare the long-time and microcanonical averages of local 


operators. The presence of chaotic initial states across 
the spectrum is discussed in Sec. VB. Results for the 
actual time evolution of local properties in the system 
are presented in Sec V C and VD. In Sec. VI, we discuss 
the relaxation time scales and the process of prethermal¬ 
ization in the system. Finally, we briefly summarize our 
findings in Sec. VH. 


II. MODEL AND QUENCH DYNAMICS 

We study thermalization in a finite isolated quantum 
system with bosons trapped in a double-well potential 
shown in Fig. 1. A single boson can occupy four modes 
(two single-particle levels in each well), which can be de¬ 
scribed by the following two-level generalization of the 
Bose-Hubbard model [32, 33]: 





FIG. 1. (Color online) Schematic of the double-well potential 
with two energy levels. The interlevel coupling [7°^ adds an¬ 
other degree of freedom in the system, which is expected to 
break integrability. 
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where bp and are the bosonic creation and annihila¬ 
tion operators of an atom in well r and energy level £, 
respectively. The Hamiltonian in Eq. (6) can be derived 
from its second-quantized form after a tight-binding ap¬ 
proximation [32]. The system is effectively mapped into 
a two-site lattice model with added dimensionality or de¬ 
gree of freedom due to the coupling with the first excited 
single-particle level. 

We consider an external harmonic potential with oscil¬ 
lator frequency wq, which is split by a focused laser beam 
located at the center of the trap and described by a Gaus¬ 
sian potential Vbexp(—a:^/2(T^). The parameters in Eq. 
(6) can be easily evaluated for a specific realization of 
the double-well potential described by the single-particle 
Hamiltonian [32], 
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The localized functions (/>f in level I G {0,1} and site 
r € are obtained by symmetric and antisymmet¬ 

ric superpositions of the eigenstates of the single particle 
Hamiltonian in Eq. (7). The corresponding eigenval¬ 
ues are = J dxcj)^*{x)Hsp4>i{x)■ The tunneling terms 
between wells are = — J dx(j)^£{x)Hsp4>^j^{x). The in¬ 
teraction terms between atoms in the same well and the 
same energy level are = g J dx\(j)f.\'^. There is also an 
interlevel coupling corresponding to the two-atom hop¬ 
ping term between energy levels in the same well given 
by = 9 J dx\(j)^{x)\‘^\4>l.{x)\'^. The limit of = 0 
is the integrable Bose-Hubbard dimer model for a dou¬ 
ble well, which was extensively studied before in different 
contexts [38-43]. Integrability of the dimer can be bro¬ 
ken by adding one or more degrees of freedom without 
changing the number of conserved quantities, which for 
the dimer model corresponds to the energy and the num¬ 
ber of bosons [44, 45]. A finite value of introduces 
additional degrees of freedom and thus breaks the inte¬ 
grability of the system. Note that upon choosing the trap 
parameters, the only free parameter in the system is the 
coupling constant g, which affects both and 17°^. This 
is a subtle detail that we need to keep in mind as blindly 
increasing the interaction g (in hopes of increasing the 
integrability breaking term will actually push the 
system towards a strongly interacting regime dominated 
by quasidegeneracies [32]. 

We now discuss the quench protocol used in our study 
of the nonequilibrium dynamics in the system. We con¬ 
sider a partitioning of the system such that each mode 
is considered as a subsystem and the remaining three 
modes will act as an environment. Initially, the num¬ 
ber of bosons in each mode is fixed; i.e., the tunneling 
and the interlevel terms are set to zero, such that the 
system starts as a product state of eigenstates of each 
subsystem. In particular, we choose a set of initial states 
corresponding to all possible Fock state configurations 
Jn) = Jrio) = The eigenstates \k) of the 

postquench Hamiltonian H with eigenvalues Ek can be 
expressed in terms of the Fock basis, 

\k) = J2c!^\n). ( 8 ) 

n 

This expansion is similar to how the mean-field or unper¬ 
turbed (in our case uncoupled) basis is physically mo¬ 
tivated in studies related to chaotic eigenstates and its 
connection to thermalization of isolated systems [28]. 

We do a sudden quench by turning on the tunneling 
and the interlevel coefficients. The wavefunction unitar- 
ily evolves in time according to 

= E Clle-^^-^l^\k). (9) 

k 

For brevity, we refer to the mean energy, which is con¬ 
served during the dynamics, 

Eo = {no\H\no)^Y.^>‘\^no\^’ (10) 

k 


as simply the initial energy of the system. 

For the postquench Hamiltonian, the barrier height is 
chosen to be Vq = 5Hu!o with width a — O.llh o, where the 
harmonic oscillator length is Iho = \/fijmuJo. The inter¬ 
action coupling g is chosen such that NU^/tux q = const. 
For this parameter space of the trap, the coefficients in 
the Hamiltonian are J^/hujo = 0.26, J^jhujo = 0.34, 
E^/HuJo = 1-25, E^/HuJo = 3.17, NU°/Hujo = const., 
= 3C/74, and = U°/2. 

We use exact diagonalization to obtain the eigenvalues 
and the eigenvectors expressed in the Fock basis. To this 
end, we obtain the DxD Hamiltonian matrix in the Fock 
basis where the size of the Hilbert space is D = {N + 
3)!/(N!)(3!) = 8436. Alternatively, the time evolution 
of the wave function can be calculated from the time- 
evolving expansion coefficients, 

IV'(i)) = (11) 

{"} 

where jn) spans all possible Fock state configurations of 
the system. Temporal evolution of the expansion coeffi¬ 
cients is computed using standard numerical integrator, 
e.g., seventh- or eighth-order Runge-Kutta method. 


III. RATIO OF CONSECUTIVE LEVEL 
SPACINGS DISTRIBUTION 

An interesting question for finite quantum systems is 
the relevance of chaoticity of the Hamiltonian in the 
study of thermalization [31]. It was conjectured using 
a system of one-dimensional lattice with spin-1/2, that 
the occurrence of thermalization depends only on the 
level of delocalization of initial states [31]. Furthermore, 
this condition has been shown to be sufficient for ini¬ 
tial states close to the middle of the spectrum irrespec¬ 
tive of integrability (or chaoticity) in both prequench and 
postquench Hamiltonian. 

The level statistics of a chaotic Hamiltonian was found 
to possess similar spectral properties as random matrices 
[46, 47]. One common measure of level statistics is the 
level spacing distribution, Sk = Ek+i — Ek, where {Ek} 
is the set of eigenenergies in ascending order [46]. In 
a Hamiltonian with time-reversal symmetry, the chaotic 
regime is identified by a level spacing distribution sim¬ 
ilar to the Gaussian orthogonal ensemble (GOE) [46]. 
The connection between nonintegrability and a GOE- 
like spectral behavior was pointed out in various works 
such as Refs. [30, 48-50]. On the other hand, integrable 
models are expected to have a Poissonian level spacing 
distribution as demonstrated in Refs. [30, 49-51]. 

Due to the large interaction strengths involved in our 
quench dynamics, i.e., NU^/> 1 for both levels i, 
it is convenient to avoid any unfolding procedure of the 
spectrum [49]. Instead, we obtain the distribution of the 
ratio of consecutive gaps between adjacent levels, Xk = 
min(sfc,Sfe_i)/max(sfc,Sfe_i) = [49, 52, 
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53]. We compare the numerically obtained distribution 
P{r) with the analytical expressions for the Poisson and 
the GOE distributions [53] 

2 2T v 

= (TTTF’ = T(l + r + r2)5/2- (12) 

A more quantitative comparison can be made using the 
mean value of r denoted as (r), which is (r)p = 21n2 —1 « 
0.3863 for the Poisson distribution and (r)GOE ~ 0.5359 
for the GOE [53]. 
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FIG. 2. (Color online) (Top) Distribution of the ratio of adja¬ 
cent level spacings r. (Bottom) Mean value (r) as a function 
of the interaction parameter. The GOE (dashed-dotted line) 
and the Poissonian (dashed line) averages are also shown. 
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There is parity or reflection symmetry about the cen¬ 
ter of the double-well system. In order to avoid mixing 
eigenenergies from different symmetry sector, we sepa¬ 
rate the eigenenergies according to its parity, i.e., even 
or odd eigenvectors in the Fock basis. We do this by 
constructing eigenstates with well-defined parity in Fock 
space [54]. For the analysis of spectral statistics pre¬ 
sented in this section, we numerically diagonalize the 
Hamiltonian in the well-defined parity basis. We denote 
the reflection operator, which exchanges the Fock state 
from left to right or vice versa as V\n^,'nPi^,n\,n}i) = 
jn^j, n° , n}j, n\). Then, we use a set of basis states given 
by: 

\n)± =-^(^\n) ±'P\n)^. (13) 

The distribution of the ratio of adjacent energy gaps 
P(r) for the even-parity sector is illustrated in Fig. 


2. First, we analyze the spectral statistics of the full 
spectrum. Notice how (r) increases with the inter¬ 
action strength before reaching its maximum around 
NU^/hujQ = 2. It then decreases as you further increase 
the interaction strength. Perhaps more interesting, the 
mean value of r becomes smaller than the Poisson distri¬ 
bution (r)p for NU^/fvjjQ > 5. This level clustering can 
be explained by the appearance of quasidegenerate eigen¬ 
vectors (in each symmetry sector) when the interaction 
terms become sufficiently larger than the tunneling 
terms [32]. Specifically, these quasidegenerate pairs 
appear as high-lying excited states of the system. Nev¬ 
ertheless, our results for the distribution of the ratio of 
consecutive level gaps imply that the postquench Hamil¬ 
tonian is nonchaotic in the chosen parameter space of 
the trap. Even though the final Hamiltonian is noninte- 
grable, the spectral statistics of the final Hamiltonian is 
actually closer to the Poisson distribution than the GOE 
for the interaction strengths considered here. This indi¬ 
cates that the system is still close to an integrable point 
after the quench. 

In the following section, we demonstrate that thermal- 
ization is still viable, in the context of the ETH, for the 
nonchaotic Hamiltonian considered here. We also argue 
later that the presence of quasidegeneracies in the spec¬ 
trum has important implications for the process of relax¬ 
ation in the system. 


IV. EIGENSTATE THERMALIZATION 
HYPOTHESIS 

The microcanonical ensemble can describe the long¬ 
time averages if the requirements of the ETH are satis¬ 
fied. Specifically, the diagonal ensemble prediction for a 

local operator (A) will be close to a value predicted by 
a thermal or microcanonical ensemble if the EEV is a 
smooth function of This requirement is fulfilled if 
the vertical width of the EEV at a chosen energy is small 
[12]. We now calculate the eigenstate expectation values 
of the mode occupation number, (fcjnf jfc), and check how 
they are distributed in the eigenenergies of the system for 
different interaction strengths. 



FIG. 3. (Color online) Distributions of EEV for NU°/huJo = 4 
but = 0. (a) hi (b) hi. 
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First, we consider the integrable limit of uncoupled 
Bose-Hubbard dimers when the interlevel coupling is ab¬ 
sent, {7°^ = 0. We plot in Fig. 3(a) the EEV of the 
dimer with the lower on-site energy while Fig. 3(b) 
shows the EEV of the dimer with higher on-site energy. 
The overall structure of the distribution of EEV appears 
similar in both dimers. That is, the expectation values 
of the mode occupation seemingly appear separated in 
flat bands, which then prohibits the EEV from being a 
smooth function of the energy. Physically, this is a direct 
consequence of NU^ / » 1. There is large range of pos¬ 

sible EEV at a given energy Ek and hence thermalization 
must be absent in this limit. 




FIG. 4. (Color online) Distributions of EEV of h% (left) and 
hi (right) for [(a),(b)] NU°/hwo = 2, [(c),(d)] NU°/tkoo = 3, 
and [panels (e),(f)] NU^/huJo = 4. Vertical lines mark 
the mean energies of initial states for the dynamics in the 
subsequent sections: (dashed) Eo/hcuo = 81.053, (dotted) 
Eo/hujo = 99.202, (solid) Eo/huo = 101.181, and (dashed- 
dotted) EojhuJo = 115.9409. 


We plot the distribution of EEV for different interac¬ 
tion strengths NU^/hujo and finite integrability breaking 
[/oi = U°/2 in Fig. 4. For NU^/hujo = 2, the dis¬ 


tribution of EEV is still regular and structured similar 
to the integrable case. Even though the EEV starts to 
clump together, separated bands are still noticeable, es¬ 
pecially in low energies. The behavior of the distribution 
of EEV, in particular its smoothness, improves as the 
interaction strength is increased from NU^/Hujq = 2 to 
NU^/Hujq = 3. For NU^/HuJo = 4, the vertical width of 
the distribution of EEV in the lower half of the spectrum 
is narrower than that in the remaining half. Note, how¬ 
ever, that eigenstates that are very close to the ground- 
state energy will always violate the ETH, at least for in¬ 
teractions considered here, since the distribution is quite 
sparse in this region. This implies the absence of ther¬ 
malization for initial energies close to the ground state. 
We point out that relaxation close to an infinite tem¬ 
perature state can be considered as a borderline case for 
which the ETH is still satisfied, as seen from the solid 
vertical line in Figs. 4(e) and 4(f). 



FIG. 5. (Golor online) Distributions of consecutive EEV gaps 
of h% for different interaction strengths. 


A more quantitative description of the smoothness of 
the EEV can be drawn from the distribution of gaps be¬ 
tween consecutive EEV [14], 


l^fel 


{k + l\d\k + 1) - {k\6\k) 


(14) 


If EEV is a smooth function of the eigenenergies, the 
distribution of jAfcj denoted by P(|Afe|) will be sharply 
peaked near jAfcj = 0. Accordingly, the width of P(|Afc|) 
decreases as the distribution of EEV smoothens. It is 
clear from Fig. 4 that the ETH works only for the bulk 
of the spectrum, i.e., in the lower half of the spectrum. 
Similar to Ref. [14], we consider only a part of the full 
spectrum by removing the lowest 10% and the highest 
25% of the eigenstates. In Fig. 5, we plot P(|Afc|) for the 
rescaled occupation number in the lower mode of the left 
well As we decrease the interaction strength, 

the distribution of jAfcj broadens and the peak decreases. 
This is consistent with the qualitative picture presented 
in Fig. 4. 

Our results in this section suggest that despite the final 
Hamiltonian being near integrability, the fact that it is 
still nonintegrable permits the fulfillment of the ETH but 
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not for all eigenstates. Therefore, the condition for an 
initial state to properly thermalize will largely depend 
on which part of the spectrum the energy of an initial 
state is located. This is indeed the case, as we see later. 


V. QUENCH DYNAMICS 

A. Comparison of diagonal to microcanonical 
expectation values 

We are interested in the relaxation dynamics if the sys¬ 
tem starts out as a product state of each mode. At long 
times, local operators in the system are expected to fluc¬ 
tuate around the mean value predicted by the diagonal 
ensemble 

jA) = J2\c!^fA,,, (15) 

k 

where = {k\A\k). Clearly, the diagonal ensemble will 
depend on the initial state through On the other 

hand, the microcanonical average depends only on the 
initial energy and for finite system will also be a func¬ 
tion of the energy window AA. The energy window is 
chosen such that it is small enough to approximate the 
mean energy of the system but still large enough to con¬ 
tain significant number of eigenstates. The appropriate 
microcanonical prediction of a local operator is obtained 
by averaging over all eigenstates within a small window 
of energy [Eq — AE, Eq + AE] 

(A)me = Jr - ^ Akk, (16) 

\Ek-Eo\<AE 

where Eg = {no\H\no) is the initial energy of the sys¬ 
tem and Ne,ae corresponds to the number of eigenstates 
within the chosen window. 



Eno / ^^0 


FIG. 6. (Color online) Relative deviations between the mi¬ 
crocanonical and the mean energy of all possible initial Fock 
states for NU°/huo = 4. 


Normally, for every initial state, one has to first calcu¬ 
late (H}me = E'ME over a range of AE and then choose 
AE that best approximates Eq. Here, in order to ease 


up on computational expense, we simply fix the energy 
window by choosing AE = 0.02Eo. Although we checked 
that the microcanonical averages obtained this way are 
still weakly dependent on other choices of AE, to keep 
track of how close the microcanonical energy is to the 
energy of the system, we define the relative deviation 

A 7^ \Eo-Eme\ 

Ah/ME.rio — ui U'J 

Aq 

The plot presented in Fig. 6, which shows the rela¬ 
tive deviations of microcanonical averages in the case 
NU^/tujjQ = 4, underpins our choice of AE. 

The postquench Hamiltonian Eq. (6) has parity or 
reflection symmetry about the center of the trap and so 
a thermalized state is expected to respect this symmetry. 
That is, the stationary properties in the left well must 
be the same as those in the right well. This is indeed 
the case for both the microcanonical and the diagonal 
ensemble predictions shown in Figs. 7 and 8. In each 
plot, both ensemble averages in the left well are sitting 
on top of the corresponding values in the right well. 
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FIG. 7. (Golor online) Diagonal ensemble (DE) vs micro- 
canonical ensemble (ME) averages of local operators for all 
possible initial Fock states, (a) n)! for NU°/hujo = 2; (b) hi 
for NU^/huJo = 2; (c) h® for NU^/hcjo = 3; and (d) hi for 
NU°/huJo = 3. 


In Figs. 7 and 8, we compare the microcanonical and 
the diagonal ensemble values calculated using Eqs. (15) 
and (16), respectively, for all initial Fock state configu¬ 
rations. Similar to the distribution of the EEV, thermal- 
ization is only viable if the distribution of the diagonal 
ensemble averages is smooth and monotonous across the 
range of possible initial energies. This is an alternative 
way of describing another signature of thermalization, 
which is the independence of the equilibrium state on 
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(c) (d) 



50 100 150 200 50 100 150 200 


FIG. 8. (Color online) Similar to Fig. 7 but NU^Ihuo = 4. 
(a) {fii)', (b) (rir); (c) ((n°)^); (d) ((n);)^). Vertical lines mark 
the energies similar to Fig. 4. 


the details of initial states. The results shown in Figs. 7 
and 8 are consistent with the prediction of thermalization 
in the context of the ETH as discussed in Sec. IV. That 
is, the microcanonical ensemble works well in predict¬ 
ing the long-time average of an initial state with energy 
somewhere in the smooth region of the EEV distribution. 
Recall that in Figs. 4(a) and 4(b), the vertical width 
of the EEV is large for NU°/Hujo = 2. This explains 
why the microcanonical averages are drastically different 
from the diagonal ensemble values for most of the initial 
states shown in Eigs. 7(a) and 7(b). Following the same 
line of reasoning, we find improvements in the behav¬ 
ior of the diagonal ensemble and its agreement with the 
microcanonical average when the interaction strength is 
increased to NU^/IujJq = 3 and NU^/Hujo = 4. These are 
shown in Figs. 7(c), 7(d), 8(a), and 8(b). We stress that 
there is a discrepancy found between the microcanoni¬ 
cal and the long-time predictions for initial states with 
sufficiently low energies, in agreement with observations 
drawn from the EEV near the ground state. 

Aside from the local one-body operator hf, we also 
calculate the ensemble averages for the local two-body 
operator ((hf)^). We plot the results in Eigs. 8(c) and 
8(d). This quantity is essential in calculating the average 
energy of a mode, which will be further investigated in 
Sec. VD. 


B. Delocalized initial states 

In the spirit of the energy shell approach, consider a 
set of unperturbed basis states \n). When a finite inter¬ 


action is turned on, this perturbation will couple some of 
these states \n). The exact eigenstates of the final Hamil¬ 
tonian can be expressed as superpositions of the unper¬ 
turbed states \k) = conjectured that 

the overlap coefficients become random variables from 
a Gaussian distribution defined by the “energy shell” 
[55, 56]. Eor systems with only two-body interaction and 
sufficiently strong perturbations, regardless of integrabil- 
ity, the shape of the energy shell follows a Gaussian pro¬ 
file centered at the initial mean energy and its width is 
equal to the energy variance [29]. Statistical description 
of the system is viable if the shape of the local density of 
states (LDoS) or strength function resembles a Gaussian 
of the same mean and variance as the energy shell [30]. 
Recently, it was conjectured that for initial states in the 
middle of the spectrum, thermalization is guaranteed if 
the initial state ergodically fills the area defined by the 
energy shell, irrespective of integrability (or chaoticity) 
of the final Hamiltonian [31]. These initial states are 
referred to as chaotic, which can be viewed as delocaliza¬ 
tion of the initial unperturbed state with respect to the 
eigenstates of the final Hamiltonian, although, whether 
this phenomenon can occur for other models and more 
importantly for initial states not too close to the center 
of the spectrum [31] is still an open question. Here, we 
respond to this query by demonstrating delocalization 
of initial states with energy away from the middle of the 
spectrum even though our postquench Hamiltonian is far 
from having a chaotic spectrum. 

The LDoS or energy distribution is defined by the dis¬ 
tribution of \C^g\'^ in the eigenvalues 

E^,iE) = J2\C^no\^S{E-E,). (18) 

k 


In practice, the LDoS is numerically obtained by divid¬ 
ing the whole spectrum of the final Hamiltonian and then 
taking in each bin [57]. In the limit of strong 

perturbation, the shape of the LDoS for finite closed sys¬ 
tem with two-body interactions follows a Gaussian dis¬ 
tribution defined by the energy shell centered at Eq and 

width a = [58-60], 

and references therein), 


W) = 


1 


\/27ra^ 


exp 


(E-Eo)^] 

2a2 


(19) 


We now investigate delocalization of initial states in 
the eigenstates of the final Hamiltonian and check if it 
coincides with the onset of thermalization found in the 
previous sections. We focus on the case NU°/hujo = 4. 
In Eigs. 4 and 8, we find smooth and monotonic behav¬ 
ior of the EEV and diagonal ensemble average in the 
lower half of the spectrum while the remaining half suf¬ 
fers from strong fluctuations. Eollowing this observa¬ 
tion, we choose four different initial energies Eq/Hujo = 
{81.05,99.20,101.18,115.94}. Specifically, we choose 
representative initial Eock states with the aforementioned 




































(a) 


C. Relaxation of mode occupation number 


(b) 





FIG. 9. (Color online) LDoS of different initial states for 
NU^/huo = 4. The dashed lines correspond to the energy 
shell. 


set of energies indicated by the vertical lines in Figs. 4 
and 8. 

We use 140 bins in the numerical calculation of LDoS 
shown in Fig. 9. We determine whether the initial state 
is chaotic with respect to the exact eigenstates by check¬ 
ing if the shape of its LDoS is close to the shape of the 
energy shell. In Figs. 9(a)-9(c), we confirm the existence 
of chaotic initial states in the system by showing that 
for these initial states, most areas defined by the energy 
shell are nicely filled. More importantly, these delocal¬ 
ized initial states display thermal behavior as shown in 
Fig. 8. Due to the system being finite, strong fluctu¬ 
ations around the mean energy is visible. Nevertheless, 
the behavior in the tails of the energy distribution is still 
quite close to the Gaussian tail of the energy shell. This is 
consistent with observations in other finite systems such 
as those found in Ref. [61]. 

Also shown in Fig. 9(d) is an initial state close to the 
middle but already in the upper half of the spectrum, 
Eq/Hujq = 115.94. Clearly, this state no longer has a 
Gaussian profile and its LDoS exhibits multiple peaks. 
It is possible that the energy distribution of this state 
can be described by a bimodal Lorentzian distribution 
similar to the one in Ref. [57] but a more careful numer¬ 
ical analysis is needed. Note that this particular case, 
Eq/Hujo = 115.94, reflects the failure of the ETH in the 
system since the distribution of the EEV for this choice 
of energy is quite broad, as seen in Eig. 4. Moreover, 
there are significantly stronger fluctuations in the diago¬ 
nal ensemble averages of initial states around this energy 
as depicted in Fig. 8. Thus, this initial state is not ex¬ 
pected to thermalize. 


Our observable of interest that is local to a subsystem 
is the occupation number of a mode hf. In this section, 
we study the relaxation dynamics from the subsequent 
time evolution of the mode occupation number (h^) after 
the proposed quench. 

In Fig. 10, we illustrate how the dynamics of the inte- 
grable dimer = 0 differ from the nonintegrable case 
with finite interlevel coupling /huo = 2. The initial 
Fock state is the same in both cases. In addition to strong 
fluctuations, recurrences occur during the time evolution 
of the mode occupation number in the integrable case. 
In contrast, we observe relaxation of the mode occupa¬ 
tion number to the value predicted by the diagonal en¬ 
semble for the nonintegrable system. Henceforth, we fo¬ 
cus on the nonintegrable case of finite interlevel coupling 
t/oi = U°/2. 



FIG. 10. (Golor online) Relaxation dynamics of occupation 
number for NU^/Twjq = 4. Integrable case U^^/huJo = 0 
(dashed line) and finite /huo = 2 (solid line). The 

diagonal ensemble average is shown in black solid lines 
for /huJo = 2. Golors denote: (red,top) n°, (ma¬ 

genta,middle) n%, (blue [gray], bottom) h\, and (green [light 
gray], bottom) h^. 

To demonstrate thermalization, or the lack thereof, we 
calculate the dynamics of (nf) for pairs of initial Fock 
states. The initial states in each pair have the same en¬ 
ergy Eq. As mentioned before, an important signature 
of thermalization is the independence of thermal predic¬ 
tions on the details of the initial state apart from con¬ 
servation laws. The conserved quantities in our system 
are the total number of bosons and the total energy. For 
our calculations, we focus on NU^/Hujo = 4 and we use 
the same set of energies as those in Sec. VB. The sys¬ 
tem is said to thermalize if the mode occupation numbers 
in each pair of initial states tend to the same stationary 
values and if these values are similar to a corresponding 
microcanonical ensemble. 

In Fig. 11(a), we present the time evolution of (nf) 
for initial states not too close to the ground-state energy. 
For this pair of initial states, the exact relaxation dynam¬ 
ics of the mode occupation numbers settle around their 
diagonal ensemble values. Moreover, the corresponding 
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500 

tUJQ 


1000 


FIG. 11. (Color online) Dynamics of occupation number for 
Fock states with the same initial energy (from top to bot¬ 
tom) Eo/hujo = 81.053, Eo/huJo = 99.202, Eo/hujo = 101.181, 
and Eo/huJo ~ 115.9409. The horizontal dashed-dotted lines 
correspond to the diagonal ensemble (DE) averages. The hor¬ 
izontal dashed lines correspond to the microcanonical ensem¬ 
ble (ME) averages. An initial product state is denoted as 


microcanonical predictions are practically indistinguish¬ 
able from the stationary value of (hf). In general, we 


find that initial states in this region of the spectrum sat¬ 
isfy initial-state independence. These observations are 
corroborated by the smooth distribution of diagonal en¬ 
semble averages in Figs. 8(a) and 8(b). 

We proceed to the characterization of the relaxation 
dynamics for initial states around the center of the spec¬ 
trum. We find that initial states in the lower half but 
closer to the center of the spectrum will still exhibit ther- 
malization. Typical results for the relaxation dynamics 
around this energy are shown in Fig. 11(b). As ex¬ 
pected for initial states in the first half of the spectrum, 
the long-time averaged occupation number in the lower 
modes are greater than those in the upper modes. This 
will continue to be valid as we increase the initial en¬ 
ergy until we finally reach the middle of the spectrum 
where the system should now relax towards an infinite 
temperature state. Due to the left-right symmetry of the 
postquench Hamiltonian, the infinite temperature state 
is characterized by having equal mean occupation num¬ 
ber between the lower modes and the upper modes, i.e., 
(fi°) = {hi) = 0.25N. This situation is quite similar to 
a two-level system, with energy spacing E 2 — Ei, cou¬ 
pled to an external heat bath at infinite temperature 
1 /T = /3 = 0 such that the ratio between the average oc¬ 
cupation numbers is ( 712 )/(ni) = exp(—/3(£'2 — Ei)) = 1. 
We present in Fig. 11(c) examples of time evolution to¬ 
wards a thermal state with almost infinite temperature. 

Even though the mode occupation numbers relax to 
their diagonal ensemble values in Fig. 11(b), there is 
a slight deviation between the microcanonical ensemble 
and the stationary values. This failure of the microcanon¬ 
ical ensemble in describing the long-time average can be 
explained by the fact that the energy distributions in the 
middle of the spectrum are quite broad, as seen in Fig. 
9(b) and Fig. 9(c). It is also possible that the discrepancy 
can be traced back from finite-size effects but further in¬ 
vestigation is required in this direction. Nonetheless, we 
loosely identify the stationary state in Fig. 11(b) as ther¬ 
mal state since the exact temporal evolutions of the mode 
occupancies still relax to the long-time average and, more 
importantly, this state respects the parity symmetry of 
the system as the bosons are evenly distributed among 
the wells. However, we point out some hints of interme¬ 
diate relaxation in Figs. 11(b) and 11(c). This is clearly 
seen from the time evolution of (h^j) in the right panel 
of the same figure. 

Examples of dynamics in the upper half of the spec¬ 
trum are presented in Eig. 11(d) for Eq/Hujo = 115.94. 
Typically, the time evolution of mode occupation num¬ 
bers will have strong temporal fluctuations in this region 
of energy. It is also quite obvious from this plot that 
the exact time evolution of the mode occupation number 
does not approach either the diagonal or the microcanon¬ 
ical ensemble averages. Instead, the steady-state values 
of the occupation number of the modes in each levels, 
are no longer equal (h^) ^ (h^). The stationary states 
in Fig. 11(d) break the reflection symmetry of the sys¬ 
tem and therefore we identify these equilibrium states as 







































10 


nonthermal. This is a concrete example of prethermal- 
ization dynamics in the system. The apparent absence 
of thermalization is elucidated by the behavior of the 
EEV shown in Fig. 4. Around Eo/fiuJo ~ 115.94, the 
distribution of (/cjnf |fc) is no longer a smooth function of 
Ek- However, it is worth mentioning that the violation of 
the ETH does not guarantee the emergence of prether- 
malized metastable states. We further investigate this 
prethermalization dynamics below. 

D. Subsystem thermalization and 
prethermalization 

In order to study subsystem thermalization, we use a 
mode partitioning scheme different from the usual spa¬ 
tial or site partitioning. For our purpose, a subsystem 
consists of only of one mode in the system. We calcu¬ 
late the dynamics of the von Neumann entropy S'vn of a 
subsystem in each level. The von Neumann entropy char¬ 
acterizes the degree of entanglement between a mode and 
the rest of the system. Moreover, this choice of partition¬ 
ing simplifies the numerical calculation of S'vn since the 
reduced density matrix of a subsystem is already diago¬ 
nal in the Fock basis. This can be seen if we write the 
many-body wave function in terms of the Fock basis 

li’it)) = ^Cn{t)\n) ( 20 ) 

n 

N 

= Y1 Cm,{p}{t)\m) ®\p), 

m—0 {p}—N—m 

where \m) is the Fock state of one mode and \p) is the 
product state of the possible occupation number in the 
other three modes. The second sum over the index {p} is 
restricted by the conservation of total number of bosons 
in the system. We can trace out the environment degrees 
of freedom \p) and obtain the reduced density matrix of 
one mode, 

H \cm,{p}it)\A\m){m\ (21) 

m^O ^{p}^N-7n ^ 

N 

= ^ \m{t)\m){m\. 

m—0 

Hence, the von Neumann entropy of a subsystem is 

N 

SvN{t) = -tr[plog/9] = - y] Xm{t)logXm{t). (22) 

m—0 

We are only concerned with the time evolution of ini¬ 
tial product state and so the von Neumann entropy will 
always start at zero. Following the quench, the von Neu¬ 
mann entropy is expected to increase as the time evolu¬ 
tion couples the modes with one another. We now study 
the dynamics of the von Neumann entropy for the same 


set of parameters and initial states used in Sec. V C. 
If the S'vn were to equilibrate, it should saturate at a 
value predicted by a Gibbs ensemble. Since the num¬ 
ber of bosons and the energy are both not conserved in 
the subsystem, the appropriate Gibbs ensemble must be 
the grand canonical one. The number of bosons and the 
energy in a subsystem are obtained using the local oper¬ 
ators 

ns=fiU Hs = U^hiini - I) -f , (23) 

respectively. Therefore, the density matrix of the grand 
canonical ensemble is 

PGC = (24) 

where Z = tr(/3Gc)- The inverse temperature /? and 
the chemical potential fi are fixed by the conditions (i) 

{Hs) = tT{HsPGc) and (ii)(ns) = tr^risPGc), where the 
overline denotes the expectation value obtained from the 

diagonal ensemble. To calculate {Hs), we need the long¬ 
time average of the two-body local operator ((nf)^) ob¬ 
tained in Sec. V A. 

Typical time evolution of the von Neumann entropy is 
shown in Fig. 12. We focus on the long-time dynam¬ 
ics and the general features of relaxation towards the 
steady-state values (or quasi-steady-state values). After 
the rapid growth in the von Neumann entropy, the behav¬ 
ior of the ensuing dynamics largely depends on the energy 
of the initial product state. For initial energy near the 
ground-state energy, the von Neumann entropy quickly 
approaches the Gibbs entropy in a single-step relaxation. 
Furthermore, the thermalization time is roughly equal to 
that of the mode occupation number. This is exemplified 
in Fig. 12(a). 

For initial states still in the lower half but closer to the 
middle of the spectrum, we observe a two-step relaxation 
process such as those presented in Figs. 12(b) and 12(c). 
The von Neumann entropy will first relax to a quasis¬ 
tationary value different from the thermal entropy. The 
metastable states manifest as prethermalization plateaus 
in the dynamics of the von Neumann entropy. We remark 
that not all of the modes will exhibit such prethermal¬ 
ization dynamics. For example, in Fig. 12(b), the modes 
in the upper levels of both wells are already fluctuat¬ 
ing close to the grand canonical entropy while the mode 
in the lower level of the right well is still trapped in a 
prethermalized state. Similar behavior is seen in Fig. 
12 (c), but here the lifetime of the prethermalized state is 
longer. The lifetime of the prethermalization plateau for 
Eo/hojo = 99.20 is two ~ 10 while for EojHuJo = 101.18 it 
is tujQ ^20. In general, we find that the time scale of the 
prethermalized regimes increases as the initial energy of 
the system increases. Gonsequently, the thermalization 
time in Fig. 12(b) {tujQ ~ 39) is smaller than that in 
Fig. 12(c) (two ~ 103). The prethermalized states in 
this part of the spectrum can be physically interpreted 
as metastable states that temporarily break the parity 
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FIG. 12. (Color online) Time evolution of the von Neu¬ 
mann entropy compared to grand canonical entropy (horizon¬ 
tal dashed lines) for (a) Eq/Twjq = 81.05, (b) Eq/TujJo = 99.20, 
(c) Eo/hujo = 101.18, and (d) Eo/hujo = 115.94. The shaded 
area emphasizes the prethermalization plateaus. The vertical 
dashed lines correspond to the thermalization time after the 
prethermalized regimes in panels (b) and (c). 


symmetry of the postquench Hamiltonian. The plateaus 
observed in the dynamics of the von Neumann entropy 
reveal the disagreement between the prethermalized en¬ 
tropy of a mode in one well and its counterpart on the 
other well. The symmetry is later restored as the system 
evolves towards parity-conserving thermal states. 

Initial states found in the upper half of the spectrum, 
such as those in Fig. 12(d), appear to get trapped in ex¬ 
tremely long-lived prethermalized states, at least within 
the numerically accessible time scales. In Fig. 12(d), the 


von Neumann entropies of each mode in the left well have 
already reached the Gibbs prediction but the modes in 
the right well still fluctuate around a nonthermal value. 
It is worth noting that the metastable states found in 
this energy region do not show clear signs of any drift 
towards thermal equilibrium. This behavior seemingly 
contradicts known results for prethermalization in other 
nearly integrable models (see Refs. [18, 26]) but there 
is no reason to rule out with certainty the possibility of 
such metastable states to decay after sufficient time. This 
issue is left for future work. 


VI. PRETHERMALIZATION AND TIME 
SCALES 

In this section, we use similar arguments as in Ref. [22] 
to understand the emergence of prethermalized states in 
the system. We split the expression for the time evolution 
of a generic local operator into three contributions, 

(i) = 5] (25) 

k^l,{k^l}^{a,b} 

+ E + E 

a^b k 

where {a, b} is a set of indices corresponding to the 
pairs of quasidegenerate eigenstates with energy differ¬ 
ence {Sab = Eb — Ea}. Written this way, it is possible 
to associate two different timescales due to the first two 
sums in Eq. (25). For brevity, we refer to each sum in 
Eq. (25), in order of its appearance, as the nondegen¬ 
erate, the quasidegenerate, and the diagonal contribu¬ 
tions. The nondegenerate contribution is associated with 
a thermalization timescale due to dephasing, as shown in 
Eq. (4). It is worth mentioning that the ETH already 
puts a strong restriction on the off-diagonal elements, 
Akk' = {k\A\k'), being small relative to its diagonal coun¬ 
terpart. Instead, we turn our attention to possible in¬ 
terplay between the structure of the spectrum and the 
coefficients . 

The energy distribution of an initial state plays an im¬ 
portant role in understanding the stages of relaxation 
dynamics in our system as they provide crucial informa¬ 
tion on the chaoticity or sparsity of the coefficients . 
The width of the energy shell measures the connectivity 
of an initial state to the eigenstates of the final Hamilto¬ 
nian. An initial state that thermalizes will never be fully 
extended with respect to the eigenstates as seen in Fig. 
9 and therefore the width of the LDoS will limit contri¬ 
butions in Eq. (25) of eigenstates with energies far from 
the initial energy. This naturally justifies a separation of 
the spectrum into different sectors. We can then analyze 
the structure of the eigenenergies within each sector. 

We divide the full spectrum into three sectors in accor¬ 
dance to the choice of initial energies in Secs. VB and 
VC; i.e., Eq/Huo = 81.05 is in the low-energy sector; 
both Eq/Hujo = {99.20, 101.18} are in the middle-energy 
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just by increasing the energy of the initial product state. 
On the other hand, most of the quasidegeneracies are 
situated in the high-energy sector. In addition to this, 
the energy difference can reach as low as ^ 10“^® in this 
sector, which is several orders of magnitude lower than 
the other sectors. These two factors combined may ex¬ 
plain why initial states in the high-energy sector exhibit 
long-lived prethermalization plateaus [see Fig. 11(d) and 
12 (d))]. 



rn 


FIG. 13. (Color online) Energy differences across the spec¬ 
trum. The eigenenergies En are arranged in increasing or¬ 
der n £ [1,13], where D is the dimension of the Hilbert 
space. (Top left) Low-energy sector n £ [1, 2000]. (Top right) 
Middle-energy sector n £ [2001,5000]. (Bottom) High-energy 
sector n £ [5001,13]. 


sector; and Eq/Hujo = 115.94 is in the high-energy sector. 
The results are presented in Fig. 13. The outliers in each 
plots represent quasidegeneracies of pairs of even and odd 
eigenstates. The energy differences, in general, decrease 
from ~ 10“® in the low-energy sector to ^ 10“® in the 
middle sector. This explains the increase in the observed 
thermalization time from t < 10^ for initial states with 
low energy to t ~ 10^ for initial energies in the midspec¬ 
trum. Also, there is perceivable increase in the prevalence 
of quasidegenerate pairs as you go to higher energies. 
This means that for sufficiently high initial energies, the 
quasidegenerate contribution in Eq. (25) becomes more 
dominant and an intermediate time scale may emerge. 
Thus, we conjecture that the presence of quasidegenera¬ 
cies led to the prethermalization observed in the system. 
Moreover, the width of the energy shell restricts the con¬ 
nectivity among the such that only the eigenenergies 
near an initial energy will contribute significantly to the 
relaxation dynamics. For this reason, initial states in the 
low-energy sector are unaffected by the quasidegeneracies 
in the high-energy sector and we only find fast and single- 
stage relaxation for such initial states [see Figs. 11(a) 
and 12(a)]. In contrast, metastable states start to ap¬ 
pear during the relatively slower relaxation dynamics of 
initial states in the middle sector [see Figs. 11(b), 11(c), 
12(b), and 12(c)]. One of the main results of this work 
is this crossover behavior from single-stage to two-stage 
relaxation process without changing the parameters of 
the Hamiltonian. Instead, we can probe this transition 


VII. SUMMARY AND CONCLUSION 

In this work, we have numerically investigated the 
relaxation dynamics following an integrability-breaking 
quench in a double-well system. We obtained the dis¬ 
tribution of consecutive level spacings in order to char¬ 
acterize the spectral statistics of the system. Then, we 
identified the postquench Hamiltonian as nonchaotic over 
a wide range of interaction parameters due to the ab¬ 
sence of level repulsion. The enhanced level clustering 
for strong interactions is attributed to the increase in the 
amount of quasidegenerate pairs. Thus, for the chosen 
trap parameters, the nonintegrable postquench Hamilto¬ 
nian is close to an integrable point. 

In order to check the requirements for the validity of 
the ETH in our system, we obtained the distribution 
of the EEV for various interaction parameters. In gen¬ 
eral, smooth distributions of the EEV are found in the 
lower half of the spectrum. By changing the interaction 
strength, we observed that the distribution of the EEV 
broadens as the number of quasidegenerate energy levels 
increases. At the moment, we only provided numerical 
evidence that the ETH holds for the system, albeit hav¬ 
ing a spectral statistics akin to integrable models. It 
would be interesting to see whether the weak ETH sce¬ 
nario studied in Refs. [16, 62] can be observed in our 
system. This can be a subject of future study involving 
careful scaling analysis of the model. 

We compared the long-time averages of local operators 
calculated using the diagonal ensemble and their micro- 
canonial values for a set of initial product states spanning 
all possible combinations of Fock states in each mode. 
The range of energy where the diagonal ensemble aver¬ 
ages are well described by the microcanonial ensemble 
is consistent with that for which the ETH is satisfied. 
Therefore, we have numerically verified that it is possi¬ 
ble for the system to thermalize. We computed the LDoS 
of typical initial states with energy within the first half 
of the energy spectrum and we classified them as chaotic 
since they ergodically fill the energy shell. These results 
allowed us to extend the scope of the main conclusion 
in Ref. [31] to include chaotic initial states away from 
the middle of the spectrum that thermalize under time 
evolution of a nonchaotic Hamiltonian. 

Using generic initial product states, we found that the 
system may exhibit thermalization in the sense that the 
exact time evolution of local operators relax to the diag- 
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onal ensemble predictions and these values are close to 
corresponding Gibbs ensemble predictions. In particu¬ 
lar, we demonstrated that the mode occupation numbers 
of certain initial states would dynamically approach the 
diagonal ensemble averages. The energy of such initial 
state is within the range of energy for which initial states 
are delocalized and the ETH is valid. Observables for ini¬ 
tial states in the high-energy region are shown to relax 
but not towards the diagonal ensemble values. This led 
us to think about the possibility of having two-stage re¬ 
laxation dynamics in the system. In fact, the emergence 
of prethermalized states became clear when we exam¬ 
ined the time evolution of the von Neumann entropy in 
each mode. Specifically, we have observed intermediate 
relaxation to metastable states for initial energies close 
to the middle of the spectrum. Finally, we have argued 
that delocalization of initial states together with the de¬ 
tails of the energy spectrum such as quasidegeneracies 
play vital roles in understanding prethermalization and 
the relaxation time scales of the system. One of the key 


contributions of this work is the possibility of probing, 
by simply changing initial product states, the continuous 
transition from single-step to two-step relaxation at fixed 
Hamiltonian parameters. 

One possible direction of future work is a careful study 
on how finite-size effects may affect the phenomenon ob¬ 
served in this work including the emergence of prether¬ 
malization and the behavior of the EEV as the size of 
the system is increased. Also, the question of whether 
an initial state with LDoS different from the energy shell 
but still has some kind of structure, as in Fig. 9(d), will 
eventually thermalize remains an open issue for future 
study. 
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